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I.  INTRODUCTION 


There  exists  a  large  class  of  vibration  problems  for  which  the  excita¬ 
tions,  or  forcing  functions,  are  random  in  nature  and,  therefore,  are  described 
statistically.  Typical  problems  include  launch  vehicle  buffeting  during  the 
transonic  time  of  flight,  random  base-shake  excitation  of  structures  or  com¬ 
ponents  during  qualification  and/or  acceptance  testing,  and  spacecraft  solar 
panel  response  analysis  for  acoustic  excitation.  The  elements  which  generally 
describe  this  class  of  problem  are: 

a.  Multi-degree-of-f reedom  dynamic  models,  resulting  in  one  or  more 
normal  modes  required  in  any  analysis. 

b.  Multiple,  random  forces  acting  on  the  system,  with  the  forces 
described  by  frequency-dependent  power  spectral  density  functions 
( PSDs ) . 

c.  Some  degree  of  correlation  between  the  forces,  with  resulting  force 
cross-correlations  that  are  nonzero. 

Traditional  analysis  procedures  (e.g.,  Refs.  1-7)  have  employed  simpli¬ 
fying  assumptions  about  one  or  more  of  the  above  problem  components.  Typical 
assumptions  have  resulted  in  analyses  which:  (1)  treat  the  system  one  mode  at 
a  time,  thus  ignoring  phasing  between  the  different  modes,  and/or  (2)  deal 
with  only  one  force  at  a  time,  thus  not  including  information  about  force 
phasing,  and/or  (3)  ignore  the  frequency-dependence  of  the  force  spectra  by 
treating  them  as  constant  throughout  certain  frequency  bands.  The  primary 
reason  for  making  these  assumptions,  it  seems,  has  been  limitations  in  compu¬ 
tational  capabilities.  With  the  advent  of  larger  and  faster  supercomputers, 
however,  the  development  of  a  procedure  for  solving  the  aforementioned  class 
of  problems,  without  making  the  traditional  simplifying  assumptions,  became 
poss  ible . 

The  development  presented  herein  is  for  linear  systems  subjected  to 
forcing  environments  made  up  of  ergodic,  random  forces,  which  are  represented 
either  as  frequency-dependent  auto-  and  cross-power  spectral  densities,  force 
time  histories,  or  a  combination  thereof.  The  analysis  procedure  is  developed 
in  a  standard  matrix  formulation,  and  is  in  essence  a  combination  of  matrix 


structural  analysis  methods  and  Gaussian  random  response  theory.  Several 
example  problems  will  illustrate  typical  applications  of  the  procedure,  as 
well  as  provide  comparisons  with  results  from  direct  numerical  solution  of  the 
equations  of  motion.  Development  of  the  procedure,  and  the  accompanying  trial 
runs  with  the  example  problems,  indicated  which  problem  parameters  need  to  be 
carefully  chosen  to  maximize  accuracy  as  well  as  computational  efficiency. 
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II.  RESPONSE  TO  RANDOM  EXCITATION  -  A  REVIEW 


The  class  of  problems  we  are  considering  involves  excitations  that  are 
random  in  nature  and,  therefore,  must  be  treated  by  statistical  means.  These 
excitations  are  nearly  ergodic  and,  therefore,  can  be  described  by  their 
frequency-dependent  auto-  and  cross-power  spectral  density  functions  (PSD  and 
cross-PSD  respectively).  This  makes  it  possible  to  solve  for  the  structural 
response  in  the  frequency  domain. 

The  behavior  of  linear  elastic  structures,  which  include  most  launch 
vehicles  and  spacecraft,  can  be  described  by  the  matrix  differential  equation 
of  motion 

[M] {x( t ) }  +  [C ] {x( t ) }  +  [K 1 {x( t ) }  =  { F ( t ) }  (1) 

In  the  above  equation  [M] ,  [C],  and  ( K ]  are,  respectively,  the  mass,  damping, 
and  stiffness  matrices,  {x(t)}  is  the  displacement  vector  (dots  denote 
differentiation  with  respect  to  time),  and  {F(t)}  is  the  vector  of 
externally  applied  forces. 

Transforming  to  normal  coordinates  we  obtain 


(I]{q(t)}  +  [2CnC0n]{q(t)}  +  [u*]{q(t)}  =  [  «4>  ]  T  { F  (  t  )  }  (2) 

where 

{x(t)}  =  [<?>]  {q( t )}  (3) 

and  [<$>]  is  the  matrix  of  mode  shape  vectors  and  (q(t)}  is  the  vector  of  modal 
displacements.  In  the  above  equations,  we  have  assumed  uncoupled  modal  damp¬ 
ing,  and  the  mode  shape  vectors  have  been  normalized  with  respect  to  the  mass 
matrix  such  that 


[<t>]T[M](4>]  =  [I] 


(M 
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By  taking  the  Fourier  transform  of  each  term  in  Eq.  (2)  and  using  the 
standard  relationships  for  Fourier  transforms  of  derivatives  we  obtain 

[  (co2  -  co-)  +  i(2^nconco)]  {q(co)}  =  {f(co)}  (5) 

In  the  above  equation,  {q(co)}  and  {f(co)}  are  the  Fourier  transforms  of  { q  ( t ) } 

T 

and  [<H  {F(t)},  respectively.  Solving  for  {q(co)},  we  obtain 


where 


{  q  (co ) }  =  [H(co)]{f(co)> 


[H(co)]  =  [(co2  -  co2)  +  i(2£  co  co)]  ^ 
n  n  n 


(6) 

(7) 


[H(co)]  is  the  structure's  modal  frequency  response  (or  admittance)  function. 
This  matrix  is  diagonal  with  a  typical  term  being 


Hk(“)  2  2 


cok  -  co  +  i2Ck“k“ 


The  physical  response  correlation  matrix  is  given  by 

r  T/2 


[R  ( x ) ]  =  lim  ^ 

x  T->oo  T 


{x(t)}{x(t  +  t)}  dt 


-T/2 


(8) 


(9) 


Substituting  (x(t)}  =  [<t>]{q(t)}  yields 

-  T/2 


[R  (t)]  =  lim  | 
x  T-*00 


[4>]{q(t)}{q(t  +  x)}T[<t>]  Tdt 


-T/2 


w 

, .  1 
lim  - 

r  T/2 

I  {q(t)}{q(t  +  x)}Tdt 

T-KO  L 

J  _x/2 

[  4>  ] 


do) 


(ii) 


=  [$i  [ Rq ( x )  ]  [<]>] 

where  [R^(t)]  is  the  modal  response  correlation  matrix. 
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It  is  shown  in  Ref.  8  that  the  auto-  and  cross-spectral  density  functions 
of  the  responses,  [S  (co)],  are  related  to  the  auto-  and  cross-spectral  density 
functions  of  the  excitations,  (Sj-(co)],  by 

[S  (co)  ]  =  [H*(u)]  (S,(oo)]  [H(w)j  (12) 

q  f 

In  Eq.  (12),  [H  (to)]  is  the  complex  conjugate  of  [  H  ( co )  ]  ,  and  [S^(co)]  and 
[  S ,-  ( co )  ]  are  the  auto-  and  cross-power  spectral  density  functions  of  the  modal 
responses  and  modal  forces,  respectively.  Since  the  Fourier  transforms  of  the 
auto-  and  cross-correlation  functions  are  the  auto-  and  cross-power  spectral 
density  functions,  respectively,  the  inverse  Fourier  transform  of  Eq.  (12)  is 
given  by 


[Rq(x)]  -  £ 


[H  (co)  ]  [Sf  (co)  ]  [H(co)  ]  eia,X  du 


To  obtain  [S^(w)]  we  begin  with  the  modal  force  correlation  matrix  given 


r  T/2 


[Rf ( c )  ]  =  lim  ^ 

T+oo  1 


<f (t  )}{f (L  +  T  )  } 1 d  t 


Substituting  {  f  ( t ) }  =  (<J>]T{F(t)}  we  obtain 


[R  (t)]  -  lim  -  I 
T+“  J 


[4>]T{F(t)}{F(t  +  x  )  }T  [  ]  d  t 


=  [ 4>  1  T [  lim  ^ 

T^o 


{F(t)}{F(t  +  t)}  dt]  [<J>] 


[*]  [Rf(t)[4»] 


where  [R  ( t  )  ]  is  the  physical  force  correlation  matrix, 
c 
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Since  [S^(u)]  is  the  Fourier  transform  of  [R^(t)],  we  obtain 


[Sf  (&>)  ] 


[<t>]T[RF(-c)]  [<|>]e  iwXdT 


=  [<t>]T[SF(a))][<|)] 


where 


[  Sp (co )  ] 


[R r(r)]e  iwxdT 
F 


(17) 


(18) 


[S  (to ) ]  is  the  physical  force  cross-PSD  matrix.  Substituting  Eq.  (17) 
r 

into  Eq.  (13)  yields 

[R  (x)]  =  <D"1([H*((o)][<t.]T[Sir(a,)][((>](H(co)])  (19) 

q  F 

where  4>  ^(  )  signifies  the  inverse  Fourier  transform  of  (  ).  Substituting 
Eq.  (19)  into  Eq.  (11)  yields 

[Rx(t)]  =  <t’~1([<t>]  [H* (<*>)]  [ 4> ] T t Sp (co ) ]  [<f>]  [H(u) ]  [(J>]T) 

Taking  the  Fourier  transform  of  both  sides,  we  obtain 

[Sx(co)]  =  [$]  [H*(u)j  [ 4> ] T [ Sp (to ) ]  [<j>]  [ H (co ) ]  [<|>]T 

This  equation  relates  the  physical  response  cross-PSD  matrix  to  the  physical 
force  cross-PSD  matrix. 


(20) 


(21) 
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The  mean  square  values  of  the  responses  are  derived  by  setting  x  =  0  in 
Eq.  (9).  Recalling  that  [R^Cx)]  and  [S^Cm)]  are  a  Fourier  transform  pair, 
and  considering  only  the  diagonal  elements  of  [S^Cu)],  we  obtain 


CO 

S^(w)  do) 

_ CO 


where  S^^(co)  signifies  the  0.-th  diagonal  element  of  [S^Cco)]. 


(22) 
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III.  NUMERICAL  IMPLEMENTATION 


For  numerical  (computer)  implementation,  we  must  deal  with  the  real  and 
imaginary  elements  of  Eq.  (21)  separately.  Thus,  we  factor  the  complex 
matrices  into  their  coincident  (real)  and  quadrature  (imaginary)  components  as 
follows : 


[Sp(cj)  ]  =  [Cp(oj)  ]  +  i  [Qp.(cj)  ] 


and 


[ H (to )  ]  =  [D(co)  ]  -  i  [E(u>)  ] 


where 


V“> 


w,  -  0) 
k 


2  2  2 
(/)Z  +  (2Ckcoku))Z 


and 


Ek<“) 


2^“k" _ 

,  2  2.2  for  ^ 

(«k  ~  «  )  +  (2^k“ku) 


(23) 


(2U) 


(25) 


(26) 


are  the  k-th  elements  along  the  diagonals  of  [D(w)]  and  [ E (to ) ] ,  respectively. 
Thus,  Eq.  (21)  can  be  written  as 


[S  (oo)]  =  [4>]{[D(w)]  +  i[E(w)]}[<j>]T{[C  (u>)3  +  i[Q  (co)]}[<j>]{[D(co)]  -  i  [E(w)  ] }  [<J>] T 
x  r  r 

(27) 


At  this  point,  for  convenience,  we  drop  the  matrix  brackets,  as  well  as  the 
(u),  and  keep  in  mind  that  all  matrices  except  [41]  are  frequency-dependent. 

Thus,  performing  the  required  algebra,  Eq.  (27)  becomes 

S  =  [ RC„R  +  RQ__,T  -  TQ_R  +  TC-T]  +  i [RQ_R  -  RC„T  +  TC„R  +  TQ„T]  (28) 

Xrrrr  t  t  t  r 


13 


J 


where 


R 


(29) 


and 


T 


(30) 


Equation  (28)  can  be  written  as 

[S  («)]  =  [C  («)]  +  i  [Q  (to)]  (31) 

XX  X 

where  [C  (w)J  and  [0  (oj ) ]  are  defined  in  Eq.  (28)  and  we  have  reintroduced  the 
x  x 

matrix  brackets. 

The  diagonal  elements  of  [ 0^(03 ) ]  are  the  auto-spectra  of  the  displace¬ 
ments.  The  cross-spectral  information  is  contained  in  the  off-diagonal  ele¬ 
ments  of  [C  («)]  and  [Q  (u)]  (the  diagonal  elements  of  [Q  (w)]  are  identi- 

XX  X 

cally  zero).  Thus,  mean  square  values  of  the  responses  of  interest  can  be 
derived  from  the  diagonal  elements  of  [C  (w)]. 

The  above  development  leads  to  displacement  response  power  spectral 
densities.  Typical  applications  can  call  for  the  derivation  of  related 
parameters  such  as  accelerations,  acceleration-based  loads,  or  combinations 
thereof.  The  power  spectral  densities  of  the  system  accelerations  can  be 
computed  by  recalling  the  relations  (Ref.  9) 

<t»(q(t))  =  iu$(q(t))  (32) 

and 

<Kq(t))  =  -  m2<J>(q(t))  (33) 
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These  alLow  us  to  write  the  acceleration-based  trequency  response  functions  as 


[ HA ( co )  ]  =  -  a)2 [ H (to )  ]  (34) 

Equations  (6)  and  (34)  relate  the  Fourier  transforms  of  the  modal  accelerations 
to  the  Fourier  transforms  of  the  modal  forces,  i.e. 

(q(to) }  =  [HA(to)l{f  (co)}  (35) 

If  power  spectral  densities  of  loads  are  required,  and  a  loads  transfor¬ 
mation  matrix,  [LTM],  which  recovers  the  loads  according  to 

{ L ( t ) }  =  [LTM]  (x(t)}  (36) 

is  available,  then  Eq.  (21)  becomes 

[SL(to)]  =  [LTM]  [<|)]  [HA*(co)]  [<t>]T[SF(to)]  [4>]  [HA(to)]  [<()]T[LTM]T  (37) 

or 

[S.  (to)]  =  [LTM]  [S--(io)  ]  [LTM]T 

-Li  X 

Another  common  application  is  the  "base-shake"  vibration  problem.  The 
extension  to  this  type  of  configuration  is  straightforward  and  is  presented  in 
Appendix  A. 

Typically,  the  number  of  forces,  N  ,  applied  to  a  system,  and  the 

r 

number  of  modes  retained  in  the  analysis,  N^,  are  significantly  less  than 
the  number  of  degrees  of  freedom  in  the  system.  Using  these  facts  when 
solving  Eq.  (21)  can  result  in  a  significant  increase  in  computational 
efficiency.  This  topic  will  be  discussed  at  length  later,  when  use  of  the 
computer  program  is  described. 

It  is  also  noted  that,  although  the  preceding  development  deals  with 
power  spectra  defined  over  both  positive  and  negative  frequency  ranges, 
practical  applications  are  limited  to  the  use  of  spectra  defined  over  a 
positive  frequency  range  only.  The  procedure  applies  directly  in  these  cases, 
as  long  as  consistency  is  maintained  throughout. 
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IV.  TEST  PROBLEMS 


To  establish  confidence  in  the  numerical  implementation  of  the  random 
response  analysis  procedure,  several  test  problems  were  solved.  The  first 
test  case  consisted  of  a  f our-degree-of-f reedom  system  (Fig.  1).  The 
structure  was  subjected  to  identical  random  forces  applied  to  the  first  two 
degrees  of  freedom.  The  force  was  derived  from  pressure  measurements  from  a 
buffet  wind  tunnel  test.  The  use  of  the  cross-spectrum  will  provide  some 
insight  into  the  importance  of  using  the  cross-spectra  in  random  response 
analyses.  Furthermore,  having  access  to  the  time  histories  of  the  forces  that 
were  used  to  produce  the  PSDs  and  cross-PSDs  allowed  for  a  direct  comparison 
of  the  subject  procedure  against  results  derived  in  the  time  domain. 

The  time  domain  analysis  results  were  obtained  by  integrating  the 
equations  of  motion  directly  with  a  fourth-order  Runga-Kutta  procedure.  The 
integration  step  was  0.001  sec,  and  the  record  length  was  30  sec.  The  mean 
square  values  of  the  responses  were  then  established  as  follows: 


2  1 
xa  “  T 


x2(t)  dt 


(38) 


where  x(t)  is  the  calculated  response  time  history,  is  the  mean 
square  value  and  T  =  30  sec.  The  frequency  domain  analysis  used  a  frequency 
increment  of  0.5  Hz,  and  the  analysis  range  was  from  0.0  to  200  Hz.  The 
frequency  increment  was  selected  to  yield  at  least  four  spectral  lines  between 
the  half-power  points  of  each  mode. 


Table  1  lists  the  root-mean-square  (RMS)  displacement  values  derived  with 
three  different  methods:  (1)  direct  numerical  integration  in  the  time  domain, 

(2)  random  response  procedure  with  cross-spectra  of  forces  included,  and 

(3)  random  response  procedure  without  cross-spectra  of  forces;  i.e.,  only  the 
auto-spectra  (diagonal  elements  of  [S  (u>)])  were  used.  The  results 

r 

indicate  that  the  responses  derived  with  the  frequency  domain  random  response 
procedure  are  in  good  agreement  with  those  from  the  time  domain  analysis,  when 
the  cross-spectrum  of  the  forces  was  included  in  the  analysis.  In  the  case 
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Figure  1.  Test  Model  Configuration 


where  only  the  auto-spectra  of  the  forces  were  used,  however,  we  see  that  the 
responses  differ  significantly — in  this  case  they  are  lower  by  up  to  27%. 

In  other  cases,  the  results  derived  by  using  only  the  force  auto-spectra  over¬ 
predicted  the  time  domain  results  (see  test  problem  2).  In  any  case,  the 
degree  of  discrepancy  between  responses  derived  using  cross-spectra  and 
responses  derived  without  using  the  cross-spectra  will  depend  on  the  degree  of 
correlation  between  the  various  forces. 

For  the  second  test  problem  the  dynamic  model  consisted  of  the  GPS 
spacecraft  coupled  to  the  Delta  II  (6925)  launch  vehicle.  Input  forcing 
functions  were  formulated  from  pressure  measurements  collected  in  a  buffet 
wind  tunnel  test  which  was  performed  using  a  rigid  scale  model  of  the  launch 
vehicle.  The  test  was  performed  at  the  Arnold  Engineering  Development  Center 
in  the  1 6T  Propulsion  Wind  Tunnel.  The  scale  model  of  the  launch  vehicle  was 
designed  and  manufactured  by  McDonnell  Douglas  Astronautics  Company  under 
contract  to  the  Air  Force  Space  Systems  Division.  In  the  test,  fluctuating 
pressure  measurements  were  taken  using  transducers  which  were  located  at 
various  stations  along  the  launch  vehicle  model.  These  transducers  were 
distributed  circumferentially  as  well  (see  Fig.  2).  The  forcing  functions 
were  derived  from  these  measured  pressure  readings  by  assigning  an  area  to 
each  transducer  (depending  on  its  location)  and  by  accounting  for  the  effect 
of  model  scale,  which  affects  both  the  amplitude  and  frequency  contents. 

* 

The  coupled  system  dynamic  model  had  47  normal  modes  to  20  Hz  ,  of  which 
3  were  rigid-body  modes  (torsional  rigid-body  motion  was  not  modeled).  The 
forces  were  applied  to  a  portion  of  the  length  of  the  vehicle  (see  Fig.  3), 
and  responses  were  recovered  at  all  degrees  of  freedom  defining  the  spacecraft 
dynamic  model.  For  illustrative  and  comparative  purposes,  we  will  concern 
ourselves  here  with  spacecraft  acceleration  results,  although  it  is  noted  that 
the  subject  procedure  was  used,  for  this  particular  system,  to  recover 
spacecraft  loads,  relative  displacements  (clearance  loss),  and  launch  vehicle 

*Delta  II  and  Titan  IV  flight  data  obtained  since  these  analyses  were  performed 
indicate  that  significant  excitation  above  20  Hz  exists.  Current  analysis 
cutoff  frequencies  are  near  40  Hz. 
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b.  TRANSDUCER  RING  LOCATIONS  ALONG  FORWARD  PORTION  OF  LAUNCH 
VEHICLE  MODEL  (indicated  by  dashed  lines) 


Figure  3.  Forcing  Function  Distribution 


accelerations  and  loads  as  well.  It  is  further  noted  that,  although  the  total 
forcing  environment  applied  to  the  coupled  system  was  formulated  from  a 
combination  of  wind-tunne l -derived  forces  (time  histories)  and  flight-data- 
based  forces  (power  spectra),  the  random  response  analysis  described  herein 
was  performed  using  the  wind-tunnel-derived  force  spectra  alone  so  as  to 
provide  a  comparison  with  results  derived  in  the  time  domain  using  the  force 
time  histories  directly. 

Figures  4a,  4b,  and  4c  show  a  typical  force  time  history,  its  power 
spectral  density,  and  one  of  its  cross-power  spectral  densities,  respectively. 
Figure  5  shows  typical  response  power  spectral  densities  at  three  locations  on 
the  spacecraft,  including  primary  structure  and  appendage  degrees  of  freedom. 
Table  2  lists  a  comparison  of  typical  resulting  RMS  values.  The  first  column 
shows  time  domain  results  obtained  with  Runge-Kutta  integration,  column  2 
shows  frequency  domain  results  obtained  with  the  random  response  analysis 
procedure  using  selected  force  cross-spectra  (i.e.,  cross-spectra  resulting 
from  forces  that  were  adjacent  to  each  other  were  included  in  the  force  PSD 
matrix,  with  the  assumption  being  that  correlation  between  forces  separated  by 
more  than  one  vehicle  station  would  have  minor  effect).  Column  3  shows  results 
derived  with  only  the  force  auto-spectra.  In  this  case  results  show  that 
omission  of  the  force  cross-spectra  will  cause  an  overprediction  of  responses 
at  some  spacecraft  locations,  and  a  slight  underprediction  at  other  locations. 


22 


IDS  /HZ 


Typical 

Cross-P 


Table  1 

Test  Problem 

Results 

R-K^ 1 ) 

RRP^2) 

RRP ( 3 ) 

*1 

0.262 

0.259 

0.213 

x2 

0.203 

0.199 

0.179 

*3 

0.177 

0.171 

0.128 

x4 

0.134 

0.131 

0.096 

Note : 

Values  are 

inches  RMS. 

(^Time  domain  Runge-Kutta  integration 
(2)prequency  domain  with  force  cross-spectra 
(3 frequency  domain  without  force  cross-spectra 


Table  2.  Typical  Spacecraft  Response  Values 


R-K^ 1 ) 

rrp<2) 

RRP ( 3 ) 

S/C  Motor  X 

Y 

0.133 

0.145 

0.134 

0.157 

0.137 

0.162 

Forward  Bulkhead  X 

Y 

0.268 

0.286 

0.276 

0.308 

0.280 

0.316 

Solar  Array  X 

Y 

0.253 

0.139 

0.258 

0.149 

0.233 

0.166 

Antenna  Tip  X 

Y 

4.65 

3.97 

4.59 

3.98 

4.20 

3.94 

Note:  Values  are  G’s  RMS 

(^Time  domain  Runge-Kutta  integration 
(^Frequency  domain  with  selected  force 
(^Frequency  domain  without  force  cross- 

cross-spectra 

-spectra 
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V.  PRACTICAL  CONSIDERATIONS 


Since  the  formulation  and  first  application  in  late  1988  of  the  procedure 
presented  herein,  there  have  been  several  Delta  II  flights  and  two  Titan  IV 
flights.  Data  obtained  from  these  flights  and  additional  review  of  the  wind 
tunnel  data  indicate  that  the  buffet  excitation  energy  content  is  still 
significant  past  the  previously  used  analysis  cutoff  frequencies  of  20  Hz. 
Studies  wich  both  vehicles  and  their  payloads  indicate  that  analysis  cutoff 
frequencies  in  the  30  Hz  to  40  Hz  range  are  needed  to  obtain  convergence  of 
spacecraft  responses.  Convergence  of  launch  vehicle  loads  is  configuration 
dependent  and  typically  occurs  below  20  Hz.  However,  convergence  studies 
should  still  be  performed  for  each  new  vehicle  configuration  and  its  payload. 

To  assure  an  adequately  conservative  prediction  of  launch  vehicle  and 
spacecraft  buffet  loads,  multiple  Mach  number  and  angle-of-attack  cases  need 
to  be  analyzed.  This  is  particularly  true  for  the  transonic  time  of  flight. 

It  is  not  unusual  to  have  different  cases  maximize  loads  in  different  areas  of 
the  launch  vehicle  and  spacecraft.  Therefore,  this  requires  that  wind  tunnel 
data  at  multiple  Mach  numbers  and  angle-of-attack  combinations  be  obtained  and 
the  appropriate  analyses  be  performed. 


An  additional  practical  consideration  is  the  analysis  frequency  incre¬ 
ment.  A  frequency  increment  that  is  too  coarse  can  result  in  a  significant 
underprediction  of  loads.  However,  it  is  also  possible  that  an 
overprediction  of  responses  can  occur.  It  has  been  found  that  adequately 
accurate  predictions  result  if  a  frequency  increment  corresponding  to  four 
spectral  lines  between  a  mode's  half-power  points  is  used.  That  is: 


where  Af  is  the  frequency  increment,  is  the  critical  damping  ratio  of 
the  mode  and  f  is  the  natural  frequency  of  the  mode. 
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As  a  finaL  note,  the  work  presented  herein  has  introduced  the  concept 
and  demonstrated  the  feasibility  of  performing  buffet  analyses  in  the  time 
domain.  This  requires  that  the  forcing  functions  be  available  in  the  form  of 
time  histories.  This  being  the  case,  rather  than  converting  the  forcing 
functions  into  the  frequency  domain,  one  may  integrate  the  equations  of  motion 
directly.  The  root-mean-square  (RMS)  values  can  then  be  computed  from  the 
response  time  histories.  Work  to  date  indicates  that  20  to  30  seconds 
of  response  data  are  needed  to  establish  RMS  values  that  adequately  charac¬ 
terize  the  statistics  of  buffet  response.  This  is  consistent  with  the  data 
collected  in  both  the  Delta  II  and  Titan  IV  wind  tunnel  tests.  However,  the 
required  lengths  of  the  time  histories  are  a  function  of  the  frequency  content 
of  the  time  histories.  Therefore,  for  each  new  configuration,  convergence 
studies  should  be  performed. 

Several  factors  make  the  calculation  of  buffet  response  in  the  time 
domain  an  attractive  proposition.  First,  the  cost  of  calculating  auto-  and 
cross-power  spectral  densities  is  avoided.  Second,  numerical  integration 
software  packages  that  directly  integrate  the  modal  equations  of  motion  are 
readily  available.  And  finally,  the  calculation  of  loads  is  simplified  since 
the  time-phasing  between  the  externally  applied  forces  and  the  acceleration  and 
displacement  responses  can  be  accounted  for  exactly  and  in  a  straightforward 
manner.  Obviously,  the  choice  of  a  time-domain  or  f requency-domain  approach 
depends  on  the  form  of  the  forcing  function  representation  and  the  respective 
costs.  However,  as  demonstrated  herein,  both  approaches  are  viable  and  each 
approach  has  its  own  advantages. 
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VI.  SUMMARY/CONCLUSIONS 


A  mul t i-f orce ,  mode-superposition  random  response  analysis  procedure  has 
been  described.  The  procedure,  which  allows  for  the  use  of  frequency-dependent 
force  power  spectral  densities  and  cross-power  spectral  densities,  is  a 
combination  of  Gaussian  random  response  theory  and  matrix  structural  analysis 
methods.  It  can  be  used  to  analyze  a  wide  variety  of  random  vibration  problems 
which  arise  in  the  structural  dynamics  area.  These  include  transonic  buffet 
loads  analyses  for  launch  vehicle/spacecraft  systems,  and  base  excitation 
analyses  typical  of  component  vibration  testing  situations.  Example  problems, 
which  help  to  establish  confidence  in  the  procedure  as  well  as  to  demonstrate 
its  use,  were  included.  Tn  addition,  the  concept  of  establishing  buffet 
responses  in  the  time  domain  was  introduced  and  its  feasibility  was  demon¬ 
strated  on  test  problems,  one  of  which  was  the  Delta  II/GPS  system. 

The  accompanying  FORTRAN  computer  program,  as  well  as  a  guide  describing 
its  use,  are  included  as  a  stand-alone  appendix.  The  program  has  been  written 
with  sufficient  flexibility  so  as  to  allow  the  user  to  adapt  it  to  a  wide 
variety  of  situations  and  needs;  however,  the  main  components  of  the  code  are 
written  with  an  eye  toward  deriving  maximum  computational  efficiency  from 
various  situations  which  often  arise  when  treating  the  aforementioned  types  of 
problems . 
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NOMENCLATURE 


[C]  =  damping  matrix 

[C  (go))  =  real  part  of  [S„(w)] 
r  r 

[C  (co)]  =  real  part  of  [S  (w)] 
x  x 

[D(co)]  =  real  part  of  [H(co)j 
[E(co)j  =  imaginary  part  of  [H(w)] 

{t(t)}  =  vector  of  modal  forces 
(F(t)}  =  vector  of  physical  forces 

[H(w)]  =  diagonal  matrix  of  modal  admittance  functions  for  displacement 

[HA(gj)]  =  diagonal  matrix  of  modal  admittance  functions  for  acceleration 

[H(o))]  =  matrix  of  physical  admittance  functions 
i  =  V - 1 

[I]  =  identity  matrix 
[K]  =  stiffness  matrix 

[M]  =  mass  matrix 
Np  =  number  of  applied  forces 

=  number  of  modes  retained  in  analysis 
{ q  f  t  ) }  =  vector  of  modal  displacements 
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{ q ( t ) }  =  vector  of  modal  velocities 

{ q ( t ) }  =  vector  of  modal  accelerations 

f  Q_(u )  ]  =  imaginary  part  of  [S_(w)] 
r  r 

[Qx(»)J  =  imaginary  part  of  [S^Cu)] 

[R  (t)]  =  correlation  matrix  of  vector  { f ( t ) } 

[R_(t)|  =  correlation  matrix  of  vector  { F ( t  ) } 

r 

[R  (x)|  =  correlation  matrix  of  vector  (q(t)} 
q 

[R  (t)]  =  correlation  matrix  of  vector  {x(t)} 

[Sp(u)l  =  cross-power  spectral  density  matrix  of  physical  forces 

[S^Cu)]  =  cross-power  spectral  density  matrix  of  physical  displacements 

x^  =  mean  square  value  of  x^ 

(x(t)}  =  vector  of  physical  displacements 

{x(t)}  =  vector  of  physical  accelerations 

(,  =  critical  damping  ratio 
2 

o  =  variance  (square  of  standard  deviation) 

[ 4>  1  -  matrix  of  normal  modes 
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r.OO, 

\y  j 


=  nnrmal  mn  H  g 


[ 4>  ]  =  niode  shape  matrix  for  force  application  degrees  of  freedom 
F 


$(  )  =  Fourier  transform  of  (  ) 


$  (  )  =  inverse  Fourier  transform  of  (  ) 


co  =  analysis  frequency  (rad/sec) 


co  =  natural  frequency  (rad/sec) 


{  }  =  complex  conjugate  of  {  } 
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APPENDIX  A 


FORMULATION  OF  EQUATIONS  OF  MOTION  FOR  BASE  EXCITATION 

For  base  excitation  where  no  relative  motion  between  the  support  points 
is  allowed,  the  equations  of  motion  for  a  linear  elastic  structure  are 

lMi<W  +  IcK*rEL>  +  [k!{xrel>  =  {0}  (A-1} 

whe  re 

=  absolute  acceleration  vector 

=  reLative  velocity 

=  relative  displacement 

and  [M|,  [C],  and  [K]  are  the  mass,  damping,  and  stiffness  matrices,  respectively 
Subs t i tut ing 

+  n-RB](xB>  (A-2) 

=  -WHV'V  (A 

where 

x  =  base  acceleration 

D 

a  i  i  f  i 

L  =  rigid  body  vectors  referenced  to  base  of  structure 
R  B 

Using  the  base- fixed  modes,  [<J>],  to  transform  to  normal  coordinates,  we  obtain 

[  I  ]  { q  >  +  [2£  co  |  { q }  +  =  —  [  r  ,  { X  }  (A-4  ) 

n  n  n  o 


^XABS^  ^XREL^ 


yields 


[MHxREL}  +  IC]{xREL 


XABS 

XREL 

XREL 
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whe  re 


UREL)  =  f  4>  1  { q}  (A-5) 

f  r]  =  C4>]t[m]  [lrb]  (  a-6  ) 

and  [<{>]  has  been  normalized  such  that 

[<$>]  T[M]  [y]  =  [I]  ( A—  7  ) 


To  enable  solution  using  the  random  response  program,  the 
vector  1  treated  as  part  of  a  set  of  mixed  coordinates,  i.e. 


base  excitation 


[ 

0 

*B 

0 

XB 

i 

,  + 

4- 

I 

q 

2£  a) 
n  n 

q 

2 

0) 

n 

q 

-r 

— 

— 

_ 

_ 

_ 

—  _ 

(A-8 ) 


Using  Eq.  (A-2),  we  can  recover  the  absolute  accelerations  in  terms  of  the 
response  quantities  from  Eq.  (A-8)  as  follows: 


— 

— 

x„ 

x„ 

fW  * 

lrb 

4> 

B 

•  =  [ATM] ' 

B 

— 

— 

q 

q 

(A-9 ) 


loads  can  now  be  recovered  with  an  acceleration  based  loads  transformation 
matrix,  [LTM],  as  follows: 


{Loads} 


[LTM]  [ATM] 


q 


(A-LO ) 


To  recover  statically-determinate  base  reaction  forces,  a  force  trans¬ 
formation  matrix,  [FTM],  can  be  used,  i.e. 

[FTM]  =  [LRg]T[M]  [ATM] 

To  assemble  the  complete  recovery  matrix,  the  separate  recovery  matrices  can 
o  .  augmented  rowwise. 


38 


APPENDIX  B 


COMPUTER  PROGRAM  USER'S  GUIDE 


1 .  BASIC  PROGRAM  DESCRIPTION 

The  random  response  procedure  (RRP)  computer  program  was  written  in 
FORTRAN  77  and  was  designed  to  run  on  CRAY  X/MP  or  CDC  mainframes.  The  code 
consists  of  a  main  program,  which  does  the  bulk  of  the  numerical  manipulation, 
and  two  subroutines — one  which  assembles  the  force  PSD  matrices  at  each 
analysis  frequency,  and  one  which  assembles  the  frequency  resoL  se  function 
matrices  at  each  analysis  frequency.  The  program  accepts  input  model  and 
forcing  spectra  data  in  a  certain  format,  and  solves  the  random  response 
equations  described  in  Sections  II  and  III  to  produce  output  in  the  form  of 
response  power  spectral  density  functions  and  root  mean  square  (RMS)  values. 

2.  INPUT  DATA 


A.  General 


The  input  data  formats  for  the  program  were  designed  to  be  consistent 
with  other  programs  and  subroutine  libraries  available  and  in  use  at  The 
Aerospace  Corporation.  These  formats,  although  possibly  in  use  only  at 
Aerospace,  are  straightforward  and  should  pose  no  problems  for  the  analyst. 


B .  Dynamic  Mode  1 

The  model  data  requirements  are  as  follows: 


Force  transformation  matrix  consisting  of  mode  shape 
vectors  transposed  and  collapsed  to  the  degrees  of 
freedom  where  forces  are  applied. 

Size  =  NMODES  x  NFORCE 
(ii)  [2CnUn]  Diagonal  modal  damping  matrix. 

Size  =  NMODES  x  NMODES 
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(iii)  [10^]  Diagonal  modal  stiffness  matrix. 

Size  =  NMODES  x  NMODES 

( iv)  [ TA j  Recovery  transformation  matrix.  Recovers  parameters  of 

interest  (loads,  physical  accelerations,  etc.)  when 
post-multiplied  by  q  or  q.  For  example,  if  accelera¬ 
tions  are  desired,  [TA]  =  [<j>]  ;  if  loads  are  desired, 
[TA]  =  [LTM]  [<t> ]  . 

Size  =  NLOAD  x  NMODES 

The  ordering  of  the  modes  should  be  from  lowest  natural  frequency  to 
highest.  Therefore,  ‘'he  rigid-body  modes  should  be  the  first  in  the  modal 
matrices.  The  matrices  should  be  on  one  permanent  file  or  tape,  specified 
locally  as  "unit  32"  in  the  Job  Control  Language  (JCL),  in  the  order  listed 
above,  and  with  end-of-file  markers  between  the  individual  matrices.  The 
matrices  must  be  written  in  that  file  in  Aerospace  Corporation  "MATRIX"  format 
(see  Appendix  C). 


C.  FORCING  FUNCTIONS 


The  forcing  functions  must  be  available  as  force  power  spectral  densi¬ 
ties  and  cross-power  spectral  densities.  These  inputs  should  be  on  a  tape  or 
permanent  file  specified  locally  as  "unit  33"  in  the  JCL.  The  default  read 
option  assumes  (for  the  case  where  cross-PSDs  are  used)  that  the  cross-PSDs 
are  written  on  the  same  file  as  the  auto-PSDs .  This  can  easily  be  altered  by 
the  user  if  he  or  she  wishes  to  use  auto-PSDs  and  cross-PSDs  that  are  written 
on  separate  files;  however,  the  frequency  spacing  of  the  two  files  must  be 
consistent.  The  frequency  band  over  which  the  PSDs  are  defined  is  arbitrary 
(but  should,  naturally,  be  kept  in  mind  when  the  analysis  bands  are  deter¬ 
mined).  The  default  force  PSD  interpolation  option  assumes  that  the  PSDs  are 
defined  on  a  linear  scale,  however,  this  can  be  changed  to  accommodate  input 
taken  directly  from  log-log  plots,  if  need  be. 

The  cross-PSDs  read  from  the  input  unit  are  mapped  to  the  appropriate 

arrays  representing  the  real  and  imaginary  components  of  the  S  matrix, 

r 

denoted  "Cp"  and  "Q^,1  in  Equation  (23)  of  the  main  text,  and  denoted  simply 
as  "C"  and  "Q"  in  the  program.  This  mapping  is  accomplished  in  the  "PSD" 


subroutine  through  the  use  of  the  arrays  IROW  and  ICOL,  which  establish  a 
correspondence  between  the  order  that  the  cross-PSD  data  is  written  on  the 
input  tape  and  the  appropriate  locations  in  the  "C"  and  "Q"  arrays. 

The  operation  of  the  "PSD"  subroutine,  at  each  analysis  frequency,  is 
as  f o 1  lows : 


i)  The  main  program  passes  the  value  of  the  analysis  frequency  to  the 
subroutine . 

ii)  The  subroutine  determines  whether  or  not  the  analysis  frequency 
lies  between  the  two  previously  defined  interpolation  limit 
frequencies.  Based  on  this  determination,  either  iiia  or  iiib  is 
executed . 

iiia)  If  the  analysis  frequency  is  within  the  limit  frequencies,  the 
input  spectra  at  those  frequencies  are  interpolated  to  the  analysis 
f  requency . 

iiib)  If  the  analysis  frequency  is  greater  than  the  previously  defined 
upper  limit,  the  next  record  of  input  spectra  is  read  from  the 
input  tape,  new  lower  and  upper  interpolation  limit  frequencies  are 
defined,  and  the  input  spectra  are  interpolated  to  the  analysis 

f  requency . 

iv)  Using  the  resulting  interpolated  input  spectra,  arrays  "C"  and  "Q" 
are  assembled  through  the  use  of  the  mapping  arrays  IROW  and  ICOL. 

v)  The  arrays  "C"  and  "Q"  are  passed  back  to  the  main  program. 


3.  OUTPUT  SPECIFICATION 


The  output  from  the  program  consists  of: 

a.  Power  spectral  density  functions  of  the  responses  defined  by  the 
rows  of  the  recovery  matrix  [TA] . 

b.  Response  RMS  values  derived  from  the  response  power  spectral 
densi t ies . 


The  output  is  written  to  a  tape  or  permanent  file  specified  locally  as  "unit 
34"  using  a  binary,  unformatted  WRITE  statement.  Thus,  for  each  analysis 
frame  defined  by  frequency  F,  the  statement 


WRITE  (34)  F,  (R(M),  M=1 ,  NLOAD ) 


is  executed.  After  the  PSD  functions  have  been  written  to  the  output  file,  an 
End-of-File  marker  is  written.  Following  the  EOF  marker,  the  RMS  values  are 
written  to  a  single  record  with  the  statement 

WRITE  (34)  (S(J),  J=l,  NLOAD) 

The  RMS  file  is  followed  by  another  EOF  marker. 

The  printed  output  consists  of: 

i)  A  summary  of  the  modal  frequencies  and  damping  values  defining  the 
input  dynamic  model. 

ii)  The  PSD  values  for  the  first  10  responses  (this  can  be  altered  by 
the  user)  at  each  analysis  frequency. 

iii)  A  summary  of  the  RMS  values  for  all  the  responses. 
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PROGRAM  NOMENCLATURE 


1 .  USER-SET  PARAMETERS 

A .  Integer 

NMODES  =  total  number  of  normal  modes  in  dynamic  model 
NRIG  =  number  of  rigid  body  modes 
NFORCE  =  number  of  applied  forces 

NCROSS  =  total  number  of  independent  force  cross-PSDs 

NFREQ  =  total  number  of  spectral  lines  in  analysis 

NLOAD  =  number  of  recovered  responses  (equal  to  the  number  of  rows  in 
[TA]) 

NPRINT  =  number  of  response  PSDs  to  be  printed  on  output  (default 
value  =  10) 

3 •  Floating  Pint 

XMUF  =  model  uncertainty  factor  (default  =  1.0) 

FZERO  =  frequency  defining  the  first  spectral  line  in  analysis 
FI  =  cutoff  frequency  of  first  band 
F2  =  cutoff  frequency  of  second  band 

F5  =  cutoff  frequency  of  fifth  band^^ 

DF1  =  frequency  increment  to  be  used  between  FZERO  and  FI 

DF2  -  frequency  increment  to  be  used  between  FI  and  F2 

DF6  -  frequency  increment  to  be  used  between  FA  and  F5^^ 

EPS  =  small  parameter  which  directs  program  to  update  frequency  incre¬ 
ment  if  analysis  frequency  is  within  EPS  of  a  band  cutoff 
frequency  (default  value  =  IE-6) 


d ^More  bands  can  be  used  simply  by  defining  the  band  cutoff  frequencies,  the 
corresponding  frequency  increments,  and  by  adding  the  appropriate  number  of 
lines  in  the  portion  of  the  program  where  the  decision  is  made  to  update  the 
frequency  increment. 


A3 


c. 


Arrays 


AMP ( NMODES ) 
IROW(NCROSS) 

ICOL(NCROSS) 


modal  amplification  factors  (default  values  =  1.0) 

cross-PSD  row  mapping  array  (used  in  conjunction  with 
ICOL(.N’CROSS)) 


cross-PSD  column  mapping  array  (used  in  conjunction 
with  I  ROW ( NCROSS  ) ) 


UU 


ARRAY  DIMENSION  GUIDE 


I.  Arrays  whose  size  is 
based  on  NMuuEa: 

AF 

7 

D 

F. 

DEN- 

TEMP 

A2ZW 

AKNMODES,  nmodes) 
A2 ( NMODES ,  NMODES) 
A3  (NMODES,  NMODES ) 
A4 (NMODES ,  NMODES) 
A5 (NMODES,  NMODES) 

II.  Arrays  whose  size  is 
based  on  NFORCE: 

SF 

SFA 

SFB 

C( NFORCE,  NFORCE) 
Q( NFORCE,  NFORCE) 

III.  Arrays  whose  size  is 
based  on  NCROSS: 

SR 

SI 

SRA 

SRB 

SIA 

SIB 

IROW 

ICOL 

IV.  Arrays  whose  size  is 
based  on  NLOAD: 

X 

S 

SA 

R 


V.  Arrays  whose  size  is  based 

on  a  combination  of  parameters 

PF(NFORCE ,  MODES ) 

TA( NLOAD,  NMODES) 


Total 

=  32 

arrays 

23 

are 

1 -d imens iona 1 

9 

are 

2-dimensional 

ANALYSIS  GUIDE 


I.  DETERMINING  ANALYSIS  BANDS  AND  FREQUENCY  INCREMENTS 

(Parameters,  NFREQ ,  FZERO ,  FI,  F2,...,  DFI  ,  DF2 , . . . ) 

Determining  the  analysis  band  cutoff  frequencies  and  associated  frequency 
increments  is  an  important  part  of  setting  up  a  random  response  analysis  run. 

As  was  previously  mentioned,  the  analyst  can  use  as  many  bands  as  he  or  she 
wishes,  or,  if  constant  frequency  steps  within  bands  are  not  required,  a 
s tead i ly-increas ing  step  size  can  be  used  from  one  end  of  the  analysis  spectrum 
to  the  other,  thus  eliminating  altogether  the  need  for  band-width  and  frequency 
increment  determination.  In  any  event,  the  parameters  that  will  dictate  the 
band  cutoff  frequencies  and  frequency  increments  will  be  the  number  of  modes 
in  the  model,  the  relative  separation  of  these  modes,  their  associated  damping, 
and  the  frequency  spectrum  over  which  the  forces  contain  energy.  In  turn,  the 
total  number  of  analysis  iterations,  NFREQ,  is  determined  by  the  value  of  the 
cutoff  frequencies  and  the  associated  increments. 

The  value  of  NFREQ  must  be  determined  by  the  analyst  after  he  or  she  has 
decided  upon  the  appropriate  cutoff  frequencies  and  increments  for  the  partic¬ 
ular  problem  at  hand.  The  primary  contributor  to  analysis  accuracy  is  the 
step  size,  in  the  neighborhood  of  a  mode,  relative  to  the  half-power  bandwidth 
of  the  mode.  The  user  must  make  certain  that  the  peaks  in  the  response  spectra 
are  adequately  characterized,  thus  allowing  their  contribution  to  the  total 
area  under  the  PSD  curve  to  be  accurately  accounted  for  when  RMS  values  are 
computed.  This  concern  is  completely  analogous  to  the  importance  of  using  an 
adequate  time  increment  when  numerically  integrating  equations  of  motion  in 
the  time  domain  using  Runge-Kutta  or  other  algorithms  (see  Ref.  10).  Of 
course,  specifying  frequency  increments  that  are  overly  fine  can  result  in 
values  of  NFREQ  that  are  inefficient,  and  can  be,  in  some  cases,  prohibitively 
expensive  from  a  computational  standpoint.  This  is  the  primary  motivation  for 
leaving  the  computation  of  NFREQ  as  a  manual  exercise  for  the  user,  as  opposed 
to  making  it  an  automatic  computation  internal  to  the  program. 


In  general,  it  has  been  found  that  an  adequate  frequency  increment  in 
the  neighborhood  of  a  mode  allows  for  four  spectral  lines  between  the  half¬ 
power  points  of  the  mode's  frequency  response  function.  Thus,  since  the 
haLf-power  bandwidth  for  a  lightly  damped  system  is  given  by 


Af,  =  2(,  f 
hp  n  n 


the  analysis  frequency  increment  should  be  no  greater  than 


Af.  =  C  f  /2 

1  n  n 


The  analyst  should  specify  a  sufficient  number  of  spectral  lines  so  that 
the  analysis  is  performed  far  enough  past  the  last  mode  of  interest  (i.e.,  the 
highest  frequency  mode  of  interest)  to  include  any  measurable  contribution 
from  the  "tail"  of  that  mode.  This  determination  will  also  depend  on  the 
mode's  frequency  and  associated  damping. 

Again,  the  requirements  for  any  particular  problem  can  vary  signifi¬ 
cantly,  and  the  above  recommendations  should  not  be  construed  as  inflexible 
constraints,  but  rather  as  general  guidelines  based  on  experience  gained  in 
applications  of  the  procedure  to  date. 

2 .  DETERMINING  FORCE  CRQSS-PSD  MAPPING  ARRAY  COEFFICIENTS 

(Parameter  NCROSS  and  arrays  IROW  and  ICOL) 

The  size  of  the  force  PSD  matrix,  ana  consequently  the  size  of  the  "C" 
and  "Q"  arrays,  is  determined  by  the  number  of  applied  forces,  NFORCE.  Thus 
the  [<$>„]  matrix  is  collapsed  to  only  those  degrees  of  freedom  at  which 
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forces  are  applied.  The  diagonal  elements  of  "C"  will  _o  •  ^  the  auto-PSDs 

(real)  of  the  forces,  whereas  the  off-diagonal  elemen  s  ...  and  contain 

the  real  and  imaginary  components,  respectively,  of  Ltie  force  cross  S'-.. 

Since  these  PSDs  and  cross-PSDs  are  typically  f req acncy-dependent ,  .  d  each 

frame,  or  spectral  line  of  data,  is  written  onto  the  input  unit  as  a  separate 
record,  a  tractable  method  for  assembling  the  "C"  and  "Q"  arrays  at  each 
analysis  frequency  is  needed.  A  scheme  involving  mapping  arrays  IROW  and  ICOL 
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(dimension  NCROSS  each)  has  been  devised-  Simply  put,  these  arrays  map  the 
cross-PSD  data  from  the  individual  records  on  the  input  unit  to  the 
appropriate  ROW/COLUMN  locations  in  arrays  "C"  and  "Q".  This  occurs  within 
the  "PSD"  subroutine,  where  the  interpolation  of  the  power  spectra  is  also 
performed . 

Presently,  the  integer  elements  of  IROW  and  ICOL  are  defined  in  a  DATA 
statement  in  the  main  program  and  are  included  in  the  common  block  CBPSD. 

This  method  of  defining  the  mapping  arrays  can  be  altered  by  the  user;  for 
instance,  the  elements  could  be  defined  by  supplemental  data,  which  would  be 
input  to  the  program  by  corresponding  READ  statements.  In  any  case,  instruc¬ 
tions  for  filling  in  the  mapping  arrays  are  more  easily  described  after 
explicitly  defining  the  format  in  which  the  auto-PSDs  and  cross-PSDs  must  be 
written  to  the  input  unit. 

Each  record  must  be  written  on  the  input  unit  using  an  unformatted 
binary  WRITE  statement,  such  as 

WRITE(33)  F,  ( 8F( K ) , K=1 , NFORCE ) ,  (SR (K ) , K=1 , NCROSS ) ,  (SI (K) , K=1 .NCROSS ) 

where 

F  -  frequency  defining  the  frame 

SF  =  array  of  force  auto-PSDs  at  F,  ordered  in  agreement  with  [<J)p] 

SR  =  array  of  real  components  of  cross-PSDs  at  F,  in  an  order 
consistent  with  IROW  &  ICOL 

SI  =  array  of  imaginary  components  of  cross-PSDs  at  F,  in  an  order 
consistent  with  IROW  &  ICOL 

If  the  analyst  is  using  cross-PSD  data  for  every  possible  force  combina¬ 
tion  (i.e.,  no  elements  of  "C"  or  "Q"  are  identically  zero)  then  the  maximum 
number  of  independent  cross-PSD  functions  will  be  needed.  This  number,  given 
by 


NCROSS  =  NFORCE (NFORCE-1 )/2 
max 


^8 


is  defined  by  the  number  of  elements  in  the  upper  triangle,  not  including  the 
diagonal  elements.  The  lower  triangle  elements  are  defined  by  the  upper 
triangle  elements  due  to  the  Hermitian  property  of  the  force  PSD  matrix. 

In  many  circumstances,  particularly  those  involving  a  large  number  of 
forces,  the  correlation,  and  thus  the  cross-PSDs,  of  forces  distant  from  each 
other,  or  in  some  other  way  unrelated,  can  be  assumed  to  be  negligible.  In 
cases  such  as  these,  the  number  of  independent  cross-PSD  functions  can  be 
significantly  less  than  NCROSS  .  Often,  in  these  cases,  the  analyst  will 
arrange  the  forces  such  that  the  resulting  force  PSD  arrays  are  banded. 

The  mapping  arrays  IROW  and  ICOL  are  needed  whenever  any  force 
cross-PSDs  are  used.  The  elements  of  the  arrays  are  defined  as  follows: 

a.  IROW  and  ICOL  are  dimensioned  as  NCROSS  (Be  sure  dimensions  are 
defined  both  in  the  main  program  and  in  the  subroutines.) 

b.  IROW(K)  and  ICOL(K)  instruct  subroutine  "PSD"  to  insert  SR(K)  and 
SI (K)  to  the  appropriate  element  locations  in  the  "C"  and  "Q"  array 
(thus  all  elements  of  IROW  and  ICOL  are  integers). 

c.  Index  K  runs  from  1  to  NCROSS. 

Example : 

i)  NFORCE  =  3 

ii)  NCROSS  =  2 

iii)  Only  cross-PSD  functions  between  forces  1  and  2,  and  2  and  3  are 
being  used 

The  resulting  mapping  arrays  would  be  defined  as  follows: 

IR0W(1)  =  1  ICOL(l)  =  2 

IR0W(2 )  =  2  ICOL( 2 )  =  3 

This  would  correspond  to  an  input  unit  written  as 

WRITE( 33 )  F,SF(1),SF(2),SF(3),SR(1),SR(2),SI(1),SI(2) 

on  a  record-by-record  basis. 
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3.  A  WORD  ABOUT  COMMON  BLOCKS  AND  THEIR  ASSOCIATED  SUBROUTINES 


The  program  uses  two  subroutines:  Subroutine  "PSD",  which  interpolates 
the  input  force  PSD  functions  to  the  defined  analysis  frequencies  and  assembles 
the  resulting  force  PSD  arrays  "C”  and  "Q",  and  subroutine  "ADMIT",  which 
assembles  the  diagonal  admittance  function  arrays  "D"  and  "E"  at  each  analysis 
frequency.  Associated  with  these  subroutines  are  common  blocks  CBPSD  and 
CBADMIT.  CBPSD  contains  variables  NFORCE,  NCROSS,  IROW,  ICOL,  and  I.  CBADMIT 
contains  variables  NMODES,  NRIG,  PI,  AF,  Z,  and  AMP. 

The  analyst  is  forewarned  to  make  certain  that 

a.  The  COMMON  statements  in  the  main  program  and  in  the  appropriate 
subroutine  are  identical. 

b.  The  arrays  used  in  common  block  memory  (IROW,  ICOL,  AF,  Z,  AMP)  are 
dimensioned  correctly  both  in  the  main  program  and  m  the 
subrout ines . 

Failure  to  heed  these  reminders  will  result  in  variable  addresses  being 
assigned  to  inappropriate  variables,  resulting  in  an  incorrect  solution. 

4 .  A  WORD  ABOUT  R IG ID  BODY  MODES ,  MODAL  AMP LIFICATION,  AND  THE  "ADMIT" 

SUBROUTINE 

Of  ten  the  dynamic  model  of  a  system  includes  one  or  more  rigid  body 
modes.  The  modal  acceleration  frequency  response  function  for  these  modes  is 
simply  the  inverse  of  their  associated  modal  masses.  However,  the  displace¬ 
ment  frequency  response  functions  will  have  a  singularity  at  u>=0,  causing 
the  computation  to  "blow-up"  at  that  frequency.  Ergo,  it  is  recommended  that 
only  the  elastic  modes  be  used  in  the  model  when  displacement  frequency 
response  functions  are  implemented. 

Modifying  the  "ADMIT"  subroutine  to  implement  displacement  frequency 
response  functions  involves  changing  the  two  lines  that  read 

D(J)  =  AMP( J )*(-( 2*P I*F)**2 )*(FJ**2-F**2 )/DEN( J) 

E(J)  =  AMP (J)*(-(2*PI*F)**2)*( 2*ZJ*F*FJ ) /DEN( J ) 
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to 

D(J)  =  AMP(J)*(FJ**2-F**2)/DEN(J) 

E ( J )  =  AMP( J )*(2*ZJ*F*FJ )/DEN( J ) 

The  analyst  might  also  choose  to  amplify  the  effect  of  certain  of  the 
modes.  This  can  be  accomplished  through  the  use  of  the  AMP  array.  AMP  simply 
provides  a  scalar  multiplier  for  each  modal  frequency  response  function.  As 
such,  the  default  setting  for  all  elements  of  AMP  is  1.0.  AMP  is  dimensioned 
to  NM0DE5,  and  the  analyst  can  alter  any  of  the  individual  elements  of  AMP  in 
the  main  program,  as  long  as  it  is  done  after  the  section  of  the  program  where 
the  arrays  are  initialized,  and  before  the  section  of  the  program  where  the 
random  response  calculation  loop  begins.  These  sections  are  clearly  denoted 
by  comment  cards. 

5.  TIMING  ESTIMATES  AND  COMPILER  CHOICE 

Execution  time  is  clearly  dependent  on  the  various  parameters  defining 
the  probLem  at  hand.  An  exact  relationship  cannot,  at  this  time,  be  estab¬ 
lished.  However,  before  running  a  full  problem,  it  is  suggested  that  the  user 
run  a  small  test  case  by  setting  NFREQ  to  a  small  number,  and  by  setting  the 
time  request  in  the  JCL  to  a  small  value.  This  will  allow  the  user  to  verify 
that  the  problem  has  been  set  up  correctly,  and  will  also  enable  him  or  her  to 
make  a  rough  estimate  of  execution  time  for  the  full  run.  These  check  runs 
can  usually  be  put  through  the  computer  relatively  quickly  due  to  the  low  time 
request . 

At  present,  the  two  FORTRAN  5  compilers  available  at  Aerospace,  for  use 
on  the  CRAY  X/MP  computer,  are  the  "CFT"  compiler  and  the  "CFT77"  compiler. 

The  CFT77  compiler  is  newer  and  has  more  ANSI  standard  features.  Both  com¬ 
pilers  have  been  used  to  run  the  random  response  code,  and  it  has  been  found 
that,  although  the  executable  code  produced  by  the  CFT77  compiler  runs  fast 
than  that  produced  by  the  CFT  compiler,  the  actual  compiling  of  the  source 
code  takes  significantly  longer  with  CFT77.  This  would  indicate  that  it  would 
be  preferable  to  use  the  CFT 7 7  compiler  for  larger  problems,  where  compile 
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time  is  a  small  portion  of  overall  CPU  time,  and  the  CFT  compiler  for  smaller 
problems,  where  compile  time  could  be  a  significant  fraction  of  the  overall 
execution  time.  Also,  if  the  executable  code  is  to  be  saved  and  used  repeat¬ 
edly,  say  with  slightly  different  inputs  in  each  case,  then  the  CFT77  compiler 
would  be  preferable.  It  is  noted  that  the  program  is  written  with  an  objective 
of  exploiting  the  vector ization  capabilities  of  the  CRAY  computer.  Both  of  the 
aforementioned  compilers  allow  for  this  vectorization  feature  to  be  jtilized. 
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APPENDIX  C 


MATRIX  TAPE  FORMAT 

General  Description 

The  standard  model  file  format  used  by  the  random  response  program  is 
termed  the  MATRIX  format.  A  MATRIX  tape  may  have  any  number  of  files,  each  in 
similar  format.  Each  file  contains  a  single  matrix  in  the  format  described 
below.  The  data  are  written  in  a  standard  default,  unformatted  form. 

File  Structure 


Assume  a  matrix  of  size  NROW  by  NCOL.  It  is  written  in  MATRIX  format  as 
f ol lows : 

Record  1  =  Header  record;  length  =  7  words 

Words  1,  2,  3,  4  =  zero 

Word  5  =  NROW 

Word  6  =  NCOL 

Word  7  =  zero 

Records  2  through  (NROW  +  1)  =  Data  records;  length  =  (NCOL  +  3) 

Word  1  =  row  number 

Word  2  -  NROW 

Word  3  =  NCOL 

Words  4  through  NCOL  +  3  -  tow  of  data 

Data  followed  by  end-of-file 

Example  file  for  matrix  of  size  (5  x  7) 

Record  10000 
Record  2137 
Record  3257 
Record  4  3  5  7 

Record  5457 
Record  6557 
EOF 


5  7  0 
Row  1  of  data 
Row  2  of  data 
Row  3  of  data 
Row  4  of  data 
Row  5  of  data 
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APPENDIX  D 

RANDOM  RESPONSE  COMPUTER  PROGRAM 
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J 


S6 


JOB , JN=D$lPT$B , T=10 , MFL=500000 . 

ACCOUNT, AC=453301 ,US=18332 . 

ACCESS, DN=FT32,PDN=*********** . 

ACCESS, DN=FT33,PDN=***** ******  . 

♦ASSIGN , DN=FT34 , RF=IW, FD=CDC , CV=ON , MBS=5120 , RS=5120 . 

CFT . 

SEGLDR,GO . 

REWIND, DN=FT34. 

♦  DISPOSE , DN=FT34 , DC=ST , DF=BB , MF=MB , T£XT= » CTASK . SAVEPF , FT34 , * 

.  ’XXXXXXXX, 10=18332, ST=PF6 . ’ . 

♦SAVE , DN=FT34 , PDN=XXXX . 

♦  EOR 

PROGRAM  RANDOM 


i-  THIS  PROGRAM  COMPUTES  RESPONSE  POWER  SPECTRAL  DENSITY  FUNCTIONS 
C  (PSD 'S)  FOR  A  LINEAR  MULTI-DEGREE-OF-FREEDOM  SYSTEM  SUBJECTED  TO 
<:  ONE  OR  MORE  FREQUENCY-DEPENDENT  FORCE  PSD’S  AND  CROSS  PSD'S. 

r 

C  A  SUGGESTED  REFERENCE  IS  ATM  NO.  89(4533-01) -55 : 

C  ’A  MULTI-MODE  RANDOM  RESPONSE  ANALYSIS  PROCEDURE’ 

C 

c  APPENDIX  ’B’  OF  THIS  REPORT  SHOULD  BE  HELPFUL  AS  AN  ANALYSIS  GUIDE 

C 

c 

o«. *«,***»«*♦**♦***♦»♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦**♦*♦**♦***♦♦♦********* 

c 

c  PETER  BROUSSINOS  AUGUST, 1989 

C 

(',„*»*****»***********♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦**♦**♦******♦********** 


c 

c 

C  THE  INPUTS  TO  THE  PROGRAM  ARE: 


1 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


1  A  DYNAMIC  MODEL  OF  THE  SYSTEM,  IN  AEROSPACE  CORPORATION 
"MATRIX"  FORMAT  (SEE  REPORT,  APPENDIX  C) .  THE  MODEL 
CONSISTS  OF 

A.  FORCE  TRANSFORMATION  MATRIX  CONSISTING  OF  MODE 
SHAPE  VECTORS  TRANSPOSED  AND  COLLAPSED  TO  THE 
DEGREES  OF  FREEDOM  WHERE  FORCES  ARE  APPLIED. 

B.  DIAGONAL  MODAL  DAMPING  MATRIX. 

C  DIAGONAL  MODAL  STIFFNESS  MATRIX. 

D  LOAD  TRANSFORMATION  MATRIX  WHICH  RECOVERS 

PARAMETERS  OF  INTEREST  (LOADS,  PHYSICAL  ACCELER¬ 
ATIONS,  ETC.)  WHEN  POST-MULTIPLIED  BY  MODAL 
DISPLACEMENTS  OR  ACCELERATIONS. 


2.  FORCE  PSD ’S  AND  CROSS-PSD’S,  IN  AEROSPACE  CORPORATION 
"TRP"  FORMAT,  WHICH  IS  DESCRIBED  IN  THE  REPORT. 


C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


THE  OUTPUT  FROM  THE  PROGRAM  CONSISTS  OF: 

1 .  POWER  SPECTRAL  DENSITY  FUNCTIONS  OF  THE  RESPONSES 
DEFINED  BY  THE  ROWS  OF  THE  RECOVERY  MATRIX. 

2.  RESPONSE  ROOT-MEAN-SQUARE  (RMS)  VALUES  DERIVED  FROM 
THE  RESPONSE  POWER  SPECTRAL  DENSITIES. 


tv*************************?******************************************* 


THIS  SECTION  CONTAINS  THE  COMMON  BLOCK  AND  ARRAY  DIMENSION  STATEMENTS 

COMMON  BLOCK  CBPSD  IS  USED  IN  THE  PSD  SUBROUTINE. 

COMMON  BLOCK  CBADM  IS  USED  IN  THE  ADMIT  SUBROUTINE. 

THE  USER  MUST  MAKE  CERTAIN  THAT  THE  LINES  ARE  IDENTICAL  IN  THE 
SUBROUTINES.  ALSO,  THE  ARRAYS  THAT  ARE  USED  IN  THE  SUBROUTINES 
MUST  HAVE  THEIR  DIMENSIONS  SPECIFIED  EXACTLY  AS  THEY  ARE  IN  THE 
MAIN  PROGRAM. 


COMMON  /CBPSD/NFORCE,NCROSS, IROW, ICOL, I 

COMMON  /CBADM/NMODES , NRIG , PI , AF , Z , AMP 

DIMENSION  AF (111) , Z (111) , D (111) , E ( 1 1 1 ) , DEN (111) , SF (22) , 

4X (1356) , S (1356) ,SA(1356) , SR  (46) , SI (46) ,SFA(22) ,SFB(22) , 

4C (22 ,22) , Q (22 , 22) , A1 (111 , 1 11) , A2 (11 1 , 111) ,R(1356) , 

4A3 (111,111) ,A4(111,111) ,PF(22,111) ,TA (1356, 111) , 

&A5 (1 1 1 , 111) , SRA (46) ,SRB(46) ,SIA(46) , 

AS IB (46) , IROW (46) , ICOL (46) ,AMP(111) ,TEMP(111) ,A2ZW(111) 

C 

C 

(;+,,,*»»•»**«»»***«***.**«*****»»***»*«***»*«»*****»**********»********* 

C 

C  THE  FOLLOWING  SECTION  IS  FOR  DEFINING  THE  FORCE  CROSS-PSD 

C  MAPPING  COEFFICIENTS.  SEE  THE  REFERENCED  REPORT  FOR  DETAILS. 

C 

C 

C 


DATA  IROW/1, 1,1, 2, 2, 3, 3, 3, 4, 4, 5, 5, 5, 6, 6, 7, 7, 7, 8, 8, 9, 9, 9, 10, 10, 
411,11,11,12,12,13,13,13,14,14,15,15,15,16,16,17,17,17,18,18,19/ 
DATA  ICOL/2, 3, 4, 3,4,4, 5,6, 5,6,6, 7, 8 ,7, 8,8,9, 10, 9, 10, 10, 11, 12, 11, 
412,12,13,14,13,14,14,15,16,15,16,16,17,18,17,18,18,19,20,19,20,20/ 


C*******************************************************************+*** 

C 

C  THE  FOLLOWING  SECTION  IS  FOR  DEFINING  THE  BASIC  PARAMETERS: 

C 
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NMUDE S=4  7 
NRIG-5 

nfgku:=?o 

N  .'ROS'  -At 
...  !  J  - !  • 

NIX*  i 
NPR1NT  =  10 
XWUL  -1.000 

C 

('****»*****»«***»*****•«•***«***»«  **»**•***»**’************************* 

c 

THIS  SL  Cl  ION  IS  FOR  DEFINING  I H£I  ANALYSIS  BAND  CUTOFF  FREQUENCIES 
AND  THE  ASSOCIATED  FREQUENCY  INCREMENTS: 


FZERJ-0 . 
F 1= .  5 
F2=5 . 
F3=10. 
F4-20. 
FS=30 


DFU.5 
DF2= .01 
DF3=  025 
DF4=0 . OSDO 
DP  5- . 10000 
DF 6= . lOOCO 


PI=3 . 141592654 
cPS= . 000001 

•  •  •  ******  *  «***«*«*:>(**«**********;***********************«*************** 
*  *  *  *  **««**jk*******w********«*********«****************>t***  +  ************ 


INHIAl  IZE  ARRAYS: 

■0  -  ■>=!  ,N;'uC.T 

\  .  r  ,  ,  . ; 

C I K ) =0 . 

N  :  K )  -  o . 

X  -0 
■  :M  ;  K ;  :-0 

'  - ',) 

AMP  ,'K  ';  0 

CONTINUE 


DO  6  K=1 , NFURCE 
SKK)=0. 

SFA (K) =0 
SFB(K) =0 . 

6  CONTINUE 

DO  1  K=1 , NCR05S 
SRA(K)=0 
SRB(K)=0 
SIA(K)=0 . 
SIB(K)=0 . 
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SR(K)=0. 

SI (K)=0 . 

7  CONTINUE 

DO  8  K=1 , NLOAD 
X (K)=0 . 

S  (K)  =0 . 

SA (K) =0 . 

R(K)=0. 

8  CONTINUE 

DO  9  K=1 .NFORCE 
DO  9  L=l, NFORCE 
C (K ,L)=0 . 
Q(K,L)=0. 

9  CONTINUE 


********************************************************************** 


THE  FOLLOWING  SECTION  ALLOWS  THE  USER  THE  OPTION  OF  SETTING  THE 
MODAL  AMPLIFICATION  COEFFICIENTS  TO  VALUES  OTHER  THAN  1.00: 


DO  850  K=N1 ,N2 
AMP (K) = 

850  CONTINUE 


PUT  MODEL  MATRICES,  LTM  IN  APPROPRIATE  ARRAYS: 


1  = 

PHIF 

....  FORCING  MATRIX 

2  = 

2*ZETA*0MEGA 

...  MODAL  DAMPING 

3  = 

(OMEGA) **2 

....  MODAL  STIFFNESS 

4  = 

TA  (OR  TC) 

....  ACCELS=TA*QDD  (L0ADS=TC*QDD) 

READ (32) 
DO  800  K= 

=1 .NMODES 

800  READ(32)  JK , JK , JK, (PF (J , K) , J=l , NFORCE) 

READ (32 , END=805) 

805  READ (32) 

DO  810  K=1 .NMODES 

READ (32 , END=810)  JK , JK , JK , (TEMP ( J) , J=1 , NMODES) 
810  A2ZW(K)=TEMP(K) 

READ (32 , END=815) 

815  READ (32) 

DO  820  K=l, NMODES 

READ (32 , END=820)  JK, JK, JK, (TEMP (J) , J=1 , NMODES) 
AF (K) =SQRT (TEMP (K) ) /  (2. PI) 

IF (K  LE.NRIG)  THEN 
Z (K)=0 . 0 
ELSE 

Z (K) =A2ZW (K)/(2*2*PI*AF(K)) 

ENDIF 

820  CONTINUE 

READ (32 , END=825) 
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n  n  n  r~,  ~  o  o  o  r~,  ~i  n  o  r~,  o  r'r'ioor-innoooo  r~>  o 


825 


READ (32) 

DO  830  K=1,NL0AD 
830  READ (32 , END=835)  JK ,  JK , JK , (TA (K , J) , J=1 , NMODES) 
835  CONTINUE 


DO  610  K=l, NMODES 
610  PRINT  7000,  K,AF(K),Z(K) 

•I********************************************************************* 

*********************************************************************** 

********************************************************+*+************ 

THIS  LOOP  DOES  THE  RANDOM  RESPONSE  CALCULATION. 

EACH  ITERATION  CORRESPONDS  TO  A  NEW  SPECTRAL  LINE. 

PARAMETER  NCASE  CAN  BE  USED  IF  ONE  WISHES  TO  RUN  MULTIPLE 
CASES  PER  RUN. 


(  DO  700  NCASE=1, 

DO  12  ITAPE=1J4 

12  REWIND  ITAPE 

DO  13  ITAPE=7 , 33 

13  REWIND  ITAPE 

NLOADl=NLOAD+l 

WRITE(l)  0,0,0,0,NFREQ,NLOAD1,0 
F=FZERO-DFl 

IF (NCROSS . EQ . 0)  GOTO  150 
DO  100  1=1 ,NFREQ 

IF  MORE  ANALYSIS  FREQUENCY  BANDS  ARE  REQUIRED,  CUTOFF  FREQUENCIES 
BEYOND  F5  AND  FREQUENCY  INCREMENTS  BEYOND  DF6  MUST  BE  DEFINED  IN 
THE  'BASIC  PARAMETERS'  SECTION,  AND  THE  APPROPRIATE  ADDITIONAL 
LINES  MUST  BE  APPENDED  TO  THE  FOLLOWING  SECTION,  OR  TO  THE 
CORRESPONDING  SECTION  IN  THE  PORTION  OF  THE  PROGRAM  FOR  NCROSS=0. 


DF=DFl 

IF (A8S (F-Fl) . LE . EPS 
IF (ABS (F-F2) . LE . EPS 
IF (ABS (F-F3) . LE.EPS 
IF (ABS (F-F4) .LE.EPS 
IF (ABS (F-F5) .LE.EPS 


.OR.  F.GT.F1)  DF=DF2 
.OR.  F.GT.F2)  DF=DF3 
.OR.  F.GT.F3)  DF=DF4 
.OR.  F.GT.F4)  DF=DF5 
.OR.  F.GT.F5)  DF=DF6 


F=F+DF 

C  PRINT  1000, I, F 

C 

DO  490  M=l, NMODES 
DO  490  K=l, NMODES 
A?(M,K)-0. 

A4 (M,K)=0. 

490  A5 (M , K) =0 . 

DO  495  M=1,NL0AD 
495  R(M)=0. 

C 
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CALL  PSD  (F,C,Q) 

CALL  ADMIT  (F,D,E) 

C 

IF (NFORCE . GT . NMODES)  THEN 
DO  500  J=1 , NMODES 
DO  500  M=l# NFORCE 

Al (M, J)=PF(M, J)  *D(J) 

500  A2 (M , J) =PF (M , J) *E(J) 

ELSE 

DO  510  M=l, NFORCE 
DO  510  J-l, NMODES 

A 1  (M ,  J)  =PF (M,  J)  *D(J) 

510  A2 (M, J) =PF (M, J) *E(J) 

ENDIF 

C 

IF (NFORCE. GT. NMODES)  THEN 
DO  520  M=1 , NMODES 
DO  520  J=l, NMODES 
DO  520  K=1 , NFORCE 
DO  520  L=l, NFORCE 

A3  (M ,  J)  =  A3  (M ,  J)  -  Al(K,M)*C(K,L)*Al(L,J) 

A4(M,J)  =  A4 (M , J)  ♦  A2 (K , M) *C (K , L) *A2 (L , J) 

520  A5 (M , J)  =  A5(M,J)  ♦  Al (K,M) *Q(K,L) *A2(L,  J) 

ELSE 

DO  525  K=l, NFORCE 
DO  525  L=l, NFORCE 
DO  525  M=l, NMODES 
DO  525  J=l, NMODES 

A3  (M ,  J)  =  A3  (M , J)  +  Al(K,M)*C(K,L)*Al(L,J) 

A4  (M , J)  =  A4 (M , J)  +  A2(K,M)*C(K,L)*A2(L,J) 

525  A5 (M ,  J)  =  A5 (M ,  J)  4  Al  (K,M)  *Q,(K ,L) *A2(L,  J) 

ENDIF 

C 

IF  (NMODES. GT.NLOADS)  THEN 
DO  530  M=1 tNLOAD 
DO  530  K=l, NMODES 
DO  530  L=l, NMODES 

530  R(M)=R(M)4TA(M,K)*(  A3 (K, L) +A4 (K ,L) 4A5(K,L) +A5 (L,K)  )*TA(M,L) 

ELSE 

DO  535  K=l, NMODES 
DO  535  L=l, NMODES 
DO  535  M=1 ,NLOAD 

535  R(M)=R(M)4TA(M,K)*(  A3 (K , L) +A4 (K , L) +A5 (K , L) +A5 (L, K)  )*TA(M,L) 

ENDIF 

C 

WRITE(l)  I ,NFREQ,NL0AD1 ,F, (R(M) ,M=1 ,NLOAD) 

WRITE(34)  F  ,  (R(M) ,M=1 ,NLOAD) 

100  CONTINUE 
GOTO  220 
C 

C**********««*********************** 

c 

C  LOOP  FOR  THE  CASE  WHERE  NO  CROSS-PSD’S  ARE  USED  (NCROSS=0) 
r 

C 

150  CONTINUE 

DO  200  1=1 ,NFREQ 
C 

DF=DF1 

IF(ABS(F-F1) .LE.EPS  .OR.  F.GT.Fl)  DF=DF2 
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IF (ARS (F - F 2) .LE.EPS  .OR.  F.GT.F2)  DF=DF3 
IF(aBS(F~F3) .LE.EPS  .OR.  F.GT.F3)  DF=DF4 
IFfABSCF-F 4}  .LE.EPS  .OR.  F.GT.F4)  DF=DF5 
ir<ABS(F-F  5)  .LE.EPS  .OR.  F.GT.F5)  DF=DF6 
C 

F=F -Dr 

C  PR IN  I  1000, 1, F 
C 

DO  200  M-l , NMODtS 
DO  2  •!;  K  .  1  ,  NMODES 
A3!M,K;---0 
290  A4 (V  ,F  j  ■ 

DO  29b' M  .1  ,NLOAD 
295  R(W)--0. 

C 

CALL.  PSD  (F,C,Q) 

CALL  ADMIT  (F,D,E) 

C 

IF  (Nr !  :0[  .GT. NMOLES)  THEN 
DO  ••  j.  1  NMODES 

DO  ,0'  V.  i  CIRCE 

'  b  .  ■ )  -  P  F  { M ,  J ) » D  ( J ) 

300  A  ,  M ,  .0 -  Pf  CM,  J)  *E  (J) 

FT  SI- 

DO  "is  M=1 , NFORCE 
DO  0 : b  J -1 , NMODES 

Al(M,  J)=PF(M,  J)*D(J) 

315  A2 (M , J) =Pl  (M , J) *E ( J) 

END  IF 
C 

IF (NFORCE. GT. NMODES)  THEN 
DO  320  M=l, NMODES 
DO  320  J=l, NMODES 
DO  320  K=l, NFORCE 

A3(M,J)  =  A3 (M , J)  +  Al(K,M)*C(K,K)*Al(K,J) 

320  A4(M,J)  =  A4(M,J)  ♦  A2 (K , M) *C (K , K) *A2 (K , J) 

ELSE 

DO  325  K=l, NFORCE 
DO  325  M=l, NMODES 
DO  325  J=l, NMODES 

A3 (M, J)  =  A3(M,J)  ♦  A1(K,M)*C(K,K)*A1(K,J) 

325  A4  (M , J)  =  A4(M,J)  ♦  A2 (K , M) *C (K , K) *A2 (K , J) 

ENDIF 

C 

IF (NMODES. GT.NLOAD)  THEN 
DO  330  M=1,NL0AD 
DO  330  K=l, NMODES 
DO  330  L=l, NMODES 

330  R(M)=R(M) +TA(M,K) * (  A3(K,L)+A4(K,L)  )*TA(M,L) 

ELSE 

DO  335  K=l, NMODES 
DO  335  L=l, NMODES 
DO  335  M=1,NL0AD 

335  R(M)=R(M)*TA(M,K)*(  A3(K ,L) +A4 (K ,L)  )*TA(M,L) 

ENDIF 
C 

WRITE(l)  I,NFREQ,NL0AD1,F, (R(M) ,M=1,N10AD) 

WRITE (34)  F,  (R(M) ,M=1 ,NLOAD) 

200  CONTINUE 


03 


00(-)0f~>(~)r)0000(~>0000r~>0 


220  CONTINUE 
ENOFILE  1 
ENDFILT  34 
REWIND  1 

f^»JHL^9f.*m*******m*******»**r**  ****<*********************************** 

«**********«************************+4******+*****+*******+*******+**** 

********************************************************************** 

THIS  SECTION  INTEGRATES  THE  RESPONSE  PSD’S  TO  PRODUCE  THE  RMS  VALUES 

NOTE:  THE  FIRST  TEN  RESPONSE  PSD’S  WILL  BE  PRINTED  AT  EACH  FREQUENCY 
STEP.  THIS  CAN  BE  CHANGED  BY  SETTING  NPRINT  TO  THE  DESIRED  VALUE, 

AND  BY  MODIFYING  THE  FORMAT  STATEMENT  (LINE  4000)  BY  CHANGING  THE 
INTEGER  MULTIPLIER  ON  THE  INSIDE  PARENTHESES. 

THE  INTEGRATION  ROUTINE  USED  IS  ’TRAPEZOIDAL’.  THIS  TOO  CAN  BE 
ALTERED  BY  CHANGING  THE  WEIGHTING  VALUES  USED.  HOWEVER,  THIS  IS  NOT 
RECOMMENDED . 


READ(l)  J1 , J2, J3 , J4 , NROW, NCOL, J5 
READ(l)  J1,J2,J3,F,  (X(J) , J=1 ,NLOAD) 

PRINT  6000 

PRINT  4000,  1 , F  , (X(J) ,J=1, NPRINT) 

DO  30  J=1 ,NLOAD 
30  S(J)=X(J) *DFl/2 . 

C30  SA(J)=X(J) *DF1/  (2  .  *  (2*PI*F) **4) 

NRW=NR0W-1 

IC=1 

DO  40  M=2 , NRW 

FPREV=F 

READ(l)  Jl , J2, J3, F,  (X(J) , J=1 , NLOAD) 

W=2 . *PI«F 

PRINT  4000,  M,F, (X(J) ,J=1, NPRINT) 

C 

IF(IC.EQ.l)  THEN 
WEIGHT=2 . 

It =2 
ELSE 

WEIGHT=2 . 

IC=1 

ENDIF 

C 

Dr-F-FPREV 

33  DO  35  J=l, NLOAD 

35  S(J)=  S(J)«-X(J) *WEICHT*DF/2 . 

C35  SA( J)=SA(J) ♦X(J) *WEIGHT*DF/ (2* (W**4) ) 

40  CONTINUE 

RFAD(l)  J1,J2,J3,F,  (X(J) , J=1,NL0AD) 

W=2 . *PI*F 
DF=F-FPREV 

PRINT  4000,  NROW,F, (X(J), J=l, NPRINT) 

DO  45  J=1 , NLOAD 
45  S(J)=  S(J)+  X(J) *DF/2 

C45  SA(J)=  SA(J) ♦  X(J) *DF/ (2* (W**4) ) 

DO  50  J=l, NLOAD 
50  S(J)=  XMUF*SQRT (S(J) ) 

C50  SA(J)=XMUF*SQRT (SA(J)) 
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r->  o  o 


PRINT  2000, NCASE , NFREQ , NDOF , NMODES , NRIG , NFORCE , NLOAD , XMUF 
DR  60  J=l, NLOAD 
PRINT  3100,  J,S(J) 

60  CONTINUE 

WRITE (34)  (S(J) ,J-1, NLOAD) 

ENDFIlE  34 
700  CONTINUE 
C 

Qtm**m*»***^****^*** **************************************************** 

r******,«, *,*«**,  +  *  +  +  +  ****  +  **m**  +  +  +  *++++***+**************************** 

Ct^tmrt^m^^tfttmmtimmf*************************************************** 

1000  FORMAT (2X, ’1=  ’,15,’  ,  F=  \Fll.5) 

2000  FORMAT (’l’,2X, ’NCASE  =  ’,13,’,  NFREQ  =  \I4,/, 

&2X/ND0F  =  ’,14,’,  NMODES  =  ’,13/ ,  NRIG  =  ’,13,’,  NFORCE  = 

*13,  * ,  NLOAD  =  ’,14,  \  MUF  =  ’,F6.2,//) 

3000  FORMAT (2X,  'RMS  OF  XDD  ’,14,’  =  ’,E12.5/  =  \F7.3,’  G’) 

3100  FORMAT (2X, ’RMS  OF  DOF  ’,14,’  =  \E12.5) 


TO  PRINT  MORE  THAN  10  RESPONSE  PSDS  ON  THE  OUTPUT,  MODIFY 
THE  FOLLOWING  FORMAT  LINE  (LINE  4000)  BY  SETTING  ’N’  TO 
C  AN  APPROPRIATE  NUMBER. 

C 

4000  FORMAT (1X,I4,F8.2, 10E12 . 4) 

C4000  FORMAT (1X,I4,F8.2,10E12.4)/,N(12X, 10E12 . 4) ) 


6000  FORMAT ( ’ 1 ’ , 3X , ’ N ’ , ’  F  ’ , 9X , ’ PSDS  OF  THE  RESPONSES  ’ ,  //) 

7000  F0RMAT(2X, 'MODAL  FREQUENCY  NO.’, 15,’  =  ’,F10.3,’  HZ,’/  DAMPING 
&RATIO  =  ’ ,F7.4) 

9000  FORMAT ('1') 

STOP 

END 

£*****i«****4r************************************************************ 
(  *********»*.ai*********************************************************** 


««****•  **************************************************************** 
(  ************************************************************* ********** 


SUBROUTINE  TO  CREATE  PSD  MATRICES  AT  EACH  F 


SUBROUTINE  PSD(F,C,Q) 

COMMON  /CBPSD /NFORCE ,NCROSS, IROW, ICOL, I 
DIMENSION  SF (22) ,SFA(22) ,SFB(22) ,SR(46) , SI (46) , 

&SRA (46) ,SRB(46) ,SIA(46) ,SIB(46) ,C(22,22) ,Q(22,22) , 

AIR0W(46) , ICOL (46) 

C 

L=NFORCE 

LC=NCROSS 

ITAP=33 

C 

IF (NCROSS . EQ . 0) GOTO  200 
C 

I F (I .GT. 1)G0T0  10 

READ(ITAP)  FA, (SFA(K) ,K=1,L) , (SRA(K) ,K=1,LC) , (SIA(K) ,K=1,LC) 
READ (ITAP)  FB , (SFB (K) , K=1 , L) , (SRB (K) , K=1 , LC) , (SIB (K) , K=1 , LC) 
10  IF(F.LT.FB)GOTO  20 
GOTO  30 
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20  DO  25  K=1 , NFORCE 

C.  *********************************** 

C  IF (SFA(K) . EQ .0 .00  .OR.SFB(K) . EQ.O.OO)THEN 

C  SF (K) =SFA (K) + ( (F-FA) * (SFB(K) -SFA (K) ) ) / (FB-FA) 

C  ELSE 

C  XXPON=LOGlO (SFA IK) ) + (  (LOGIO(F) -LOGIO(FA)  )* 

C  &  (  LOGlO(SFB(K)  ) -LOG10(SFA(K)  )))/(  LOGIO(FB) -LOGIO(FA)  ) 

C  SF (K) =10 . **XXPON 

C  ENDIF 

C  25  CONTINUE 

C********«************************** 

c 

25  SF (K) =SFA (K)  +  ( (F-FA) *  (SFB (K) -SFA (K) ) ) / (FB-FA) 

DO  27  K=1 ,  LC 

SR (K) =SRA (K)  ♦  ( (F-FA) * (SRB (K) -SRA (K) ) ) / (FB-FA) 

27  SI (K) =SIA (K) ♦ ( (F-FA) * (SIB(K) -SIA (K) ) ) / (FB-FA) 

GOTO  100 
30  FA=FB 

DO  35  K=l, NFORCE 
35  SFA(K)=SFB(K) 

DO  37  K=1,LC 
SRA(K)=SRB(K) 

37  SIA(K)=SIB(K) 

READ (ITAP)  FB , (SFB (K) , K=1 , L) , (SRB (K) , K=1 , LC) , (SIB (K) , K=1 , LC) 
GOTO  10 
100  CONTINUE 

C 

DO  40  K=l, NFORCE 
40  C  (K,K)=SF  (K) 

C 

DO  45  K=1 , LC 
C ( IRQ A  (K) , ICOL (K) ) =SR (K) 

Q (IROW(K) , ICOL (K) ) =SI (K) 

C (ICOL (K) , IROW(K))=SR(K) 

45  Q (ICOL (K) , I ROW (K) ) =-SI (K) 

GOTO  099 


C  THIS  SECTION  IS  FOR  THE  CASE  WHERE  NO  FORCE  CROSS-PSDS  ARE 
0  BEING  USED  (NCROSS  =  0) . 

200  CONTINUE 

IF (I . GT . 1) GOTO  510 
READ (ITAP)  FA,(SFA(K),K=1,L) 

READ (ITAP)  FB, (SFB(K) ,K=1,L) 

510  IF(F.LT.FB)GOTO  520 
GOTO  530 

520  DO  525  K=1 , NFORCE 

C*******»*************************** 

C  IF (SFA (K) .EQ. 0.00  . OR . SFB (K) . EQ . 0 . 00) THEN 

C  SF (K)=SFA (K) ♦ ( (F-FA) * (SFB (K) -SFA (K) ) ) / (FB-FA) 

C  ELSE 

C  XXPON=LOGlO(SFA(K)) + (  (LOGIO(F) -LOGIO(FA)  )* 

C  4  (  LOGlO(SFB(K)  ) -LOG10(SFA(K)  )))/(  LOGIO(FB) -LOGIO(FA)  ) 

C  SF (K) =iO . **XXPON 

C  ENDIF 

C525  CONTINUE 

r**********  ************************* 

r 

525  SF (K)=SFA(K)  -*-((F-FA)  *  (SFB (K) -SFA (K)))  /(FB-FA) 
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GOTO  6uO 
520  FA-FB 

DO  535  K=1 , NFORCL 
535  SFA(K)=SFB (K) 

READ (I TAP)  FB, (SFB(K) ,K=1,L) 
GOTO  510 
600  CONTINUE 
C 

DO  540  K=1 , NFORCE 
540  C(K,K)=SF(K) 

C 

999  RETURN 
END 


+  +  *  +  +  +  +  *  +  4  +  +  +  +  *  +  +  +  +  +  +  +  +  +  +  +  +  + 

+  +  +  +  +  *  +  +  +  +  +  *  +  m  +  *  +  +  m  +  +  +  m  + 

SUBROUTINE  TO  CREATE  ADMITTANCE  MATRICES  AT  EACH  F 

SUBROUTINE  ADMIT (F,D,E) 

COMMON  /CBADM/NMODES ,NRIG,PI,AF,Z, AMP 

DIMENSION  DEN(lll)  ,0(111)  ,E(111)  ,AF(111)  ,Z(1U)  , AMP (111) 

C 

DO  100  J=1 , NMODES 
FJ=AF ( J) 

ZJ=Z(J) 

IF(J.LE.NRIG)  GOTO  90 
C  PRINT  999,  J , FJ , Z J 

C  999  FORMAT  (IX ,  ’  J=  ’,14/  FJ=  ’  ,F7 .3,  ’  ZJ=  \F6.3) 

DEN ( J) =((FJ**2-F**2)*+2+(2*ZJ*FJ*F)**2)*(2*PI)*+2 
IF (J . EQ . 9) DEN ( J) =0 . 

IF  (DEN ( J)  EQ.O.)  THEN 
110  PRINT  200, J,DEN(J) 

STOP 

ENDIF 

C  *********************************** 

r 

C  THH  FULL  OWING  2  LINES  SHOULD  BE  ADJUSTED  ACCORDING  TO  THE  REPORT 

IF  THt  USER  WISHES  TO  USE  DISPLACEMENT  FREQUENCY  RESPONSE  FUNCTIONS 
INSTEAD  OF  ACCELERATION  FREQUENCY  RESPONSE  FUNCTIONS. 


D(J) -AMP(J) « (- (2*PI+F) *  +  2) * (FJ++2-F**2) /DEN ( J) 
E ( jj  =AMP ( J) »(-(2*PI*F)**2)* (2*ZJ*F*FJ) /DEN(J) 


GOTO  100 
90  D(J)=1  . 

E(J)=0. 

100  CONTINUE 
125  CONTINUE 

200  F0RMAT(2X, »J=  '15, 2X, 'DEN(J)=  ’F10.5, ’  THE  DENOMINATOR  IN  THE 
^FREQUENCY  RESPONSE  FUNCTION  =  0’ ,//,30X,  ’CHECK  BASIC  PARAMETERS’, 
&/// ,23X , ’*******  EXECUTION  TERMINATED  ****.**’,/////) 

RETURN 

END 

*EOR 
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